On the local nature and scaling of chaos in weakly nonlinear disordered chains 
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The dynamics of a disordered nonlinear chain can be either regular or chaotic with a certain 
probability. The chaotic behavior is often associated with the destruction of Anderson localization 
by the nonlinearity. In the present work, it is argued that at weak nonlinearity chaos is nucleated 
locally on rare resonant segments of the chain. Based on this picture, the probability of chaos 
is evaluated analytically. The same probability is also evaluated by direct numerical sampling of 
disorder realizations, and quantitative agreement between the two results is found. 



I. INTRODUCTION 

Dynamics of classical disordered nonlinear chains is 
governed by an interplay of two fundamental phenomena: 
Anderson localization (AL) and chaos. AL, originally 
introduced as full suppression of diffusion by interfer- 
ence in random lattices [l[ and used to explain electronic 
metal-insulator transitions, turned out to be a univer- 
sal wave phenomenon and was observed in such diverse 
systems as microwaves, light, acoustic waves, and ultra- 
cold atoms @. It is most pronounced in one dimension, 
where all eigenmodes of a disordered linear system _are 
exponentially localized by arbitrarily weak disorder 

By now, AL is well understood for linear systems 
In the presence of a nonlinearity, the situation is more 
complicated, and many fundamental questions are still 
open. For example, it is unclear what happens to an ini- 
tially localized wave packet at very long times (see Ref. [|| 
for a review). In a linear system with AL, it remains ex- 
ponentially localized at all times. With nonlinearity, the 
wave packet width was found to increase as a subdiffusivc 
power of time t in numerical simulations In con- 

trast, a rigorous argument shows that at long times the 
spreading, if any, must be slower than any power of time 
10]. Analysis of perturbation theory in the nonlinearity 
suggests that there is a front propagating as In t beyond 
which the wave packet is localized exponentially [ll| . An 
indication for a slowing down of the power law has also 
been seen in the scaling analysis of numerical results [l2j . 
Finally, a possible mechanism of breakdown of subdiffu- 
sion at long times is presented in Ref. 



13] 



One of the difficulties in describing nonlinear system is 
that their dynamics can be chaotic |14l - [l6j . For Hamilto- 
nian systems, close to integrable and with a finite num- 
ber of the degrees of freedom, the volume of the chaotic 
phase space is small as long as the integrability-breaking 
perturbation (in our case, the coupling between the os- 
cillators or the nonlinearity) is weak, as guaranteed by 
the Kolmogorov-Arnol'd-Moser (KAM) theorem. 

For a random system, the description of chaotic and 
regular motion has to be probabilistic. In the pioneering 
work [l7[, the existence of a dense set of regular trajecto- 
ries was proven for a class of disordered weakly-nonlinear 
lattices. It has been argued for the disordered nonlinear 
Schrodinger chain (DNLS) that the measure of such set 
in the phase space is finite jl8]. This measure, aver- 



aged over the disorder, was recently estimated for DNLS 
both analytically [l9| and numerically (2(j. In particu- 
lar, when the nonlinearity is weak, chaos should appear 
locally. Namely, for a sufficiently long chain, arbitrarily 
separated in two segments, the probability to be on a reg- 
ular trajectory is given by the product of the two individ- 
ual probabilities for the segments. This naturally follows 
from two conditions: (i) if any of the segments is chaotic, 
the whole chain is chaotic, and (ii) the probabilities for 
the two segments are independent. Moreover, in Ref. (l9j 
an explicit mechanism for chaos generation was proposed: 
two coinciding resonances for rare combinations of three 
oscillators. Still, the results of Refs. 0, did not 
agree quantitatively, and the origin of the discrepancy 
is currently not understood. Evidence for a local origin 
of chaos has also been found in the simulations of the 
dynamics of a classical spin chain (2l| . 



In the present work, the probability of chaotic behavior 
(appearance of a non-zero Lyapunov exponent) is studied 
for another system of coupled nonlinear oscillators [see 
Eq. ([T]) below] , which turns out to be simpler to analyze 
than the DNLS. The probability of chaos is calculated at 
low energy densities e by two different methods: (i) by the 
analysis of the phase space of an effective Hamiltonian 
describing a resonant triple of oscillators, in the spirit of 
Ref. [l9j , and (ii) by direct numerical sampling of many 
disorder realizations, and counting those with nonzero 
Lyapunov exponent, analogously to Ref. 20]. 



Our main results are the following, (i) The two cal- 
culations agree quantitatively, including both the lead- 
ing low-e scaling exponent, and the numerical pref actor, 
thereby confirming the dominant role of resonant triples 
in generation of chaos at low energies, (ii) This agree- 
ment sets in at unexpectedly small values of e, while at 
moderately low e the numerical results fall on an interme- 
diate asymptotics. This intermediate asymptotics is not 
controlled by any small parameter and seems to exist for 
purely numerical reasons; the microscopic mechanism re- 
sponsible for it remains unclear at the moment. 
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II. STATEMENT OF THE PROBLEM 

Here we study the model denned by the Hamiltonian 

H{{p n ,q n }) = H ({p n ,q n }) + H int ({p n ,q n }) = 



^\2m 2 



where q ni p n are the coordinate and momentum of the 
nth oscillator, and ui 2 are independent random variables 
uniformly distributed in the interval 



W 2 



< OJZ < 



3W 2 



(2) 



W being the disorder strength. This model belongs to 
the class of models considered in Ref. [ItJ ■ It has also 
been studied in Refs. [HI, however, in the latter two 
works focused on spreading of an initially localized wave 
packet, rather than on the probability of chaos. 

The model of Eq. ([1} is different from DNLS in two 
aspects. First, it has only one conserved quantity (en- 
ergy) in contrast to the two (energy and total norm) for 
DNLS. Second, it has only one dimensionless parameter, 



Pi 



i 2 W 4 \2m 



9 



(3) 



controlling both the nonlinearity and the coupling be- 
tween the oscillators. In DNLS they are controlled sepa- 
rately by two independent dimensionless parameters, and 
when both are small, it matters which one is smaller. In 
the FSW case, both limits of weak coupling and weak 
nonlinearity correspond to e n — > 0, which makes it easier 
to analyze than DNLS. From now on, we me asure mo- 
menta, coordinates, and time in the units of \J m 3 W 4 / g, 
\J mW 2 jg, and 1/W, respectively, which is equivalent to 
setting m, W, g = 1. 

For a given realization of disorder {w 2 } and a given 
initial condition {p ni <?«}, the system trajectory is cither 
regular (quasi-periodic) or chaotic. Characterizing each 
initial condition by the typical energy e per oscillator, 
we define the probability for the initial condition to be 
chaotic as 



P(e,L)= I 9({p n ,g„},K 2 l })x 



pI 



dp n dq n 
2tc/u)„ 



(4) 



where Q({p n ,q n },{uj 2 }) = 1 if the trajectory is chaotic 
and zero otherwise. Eq. (QJ corresponds to fixing the 
energies of all oscillators to be e, and will be referred 
to as the fixed-e ensemble. Alternatively, one can fix 
the energies only on average, replacing the (5-function in 
Eq. (0} by (l/e)exp[-(p£ + ui^ql) / (2e)], i. e. a thermal 
distribution with temperature e (provided that e <C 1, 
which is our main focus) (23j. In the present paper we 



will work with the fixed-e ensemble which is easier to 
handle numerically. The corresponding initial conditions 
can be represented as 



Pn 



2e cos < 



2e 



(5) 



where the phases <p n are uniformly distributed on the 
interval [0, 27r]. 

The property of locality, mentioned in the introduc- 
tion, leads to the dependence P(e, L) = 1 — e^ w ^ L at 
sufficiently large L, where the quantity w(e) , which we 
call average chaotic fraction (as in Ref. 19]), does not 
depend on L. The main goal of the present work is to 
establish its asymptotic behavior w(e — > 0), by two meth- 
ods: (i) relating it to the chaotic phase volume of three 
resonant oscillators, as in Ref. [l9j], and (ii) by direct nu- 
merical sampling, analogously to Ref. [20]. The latter 
method also provides a check for the locality hypothesis: 
if indeed P(e, L) = 1 - e~ w ^ L , the quantity 



1 



P(e,L) 



does not depend on L, so that all curves for different L 
should collapse on a single curve w(e). Thanks to the 
relative simplicity of the system defined by Eq. (flj , both 
calculations can be carried out all the way to the final 
result, which turns out to be 

w(e^0)=Ae 2 , (6) 

where A as 1.37 • 10 3 for the fixed-e ensemble. 



III. CALCULATION 

The reduction of the problem to a few-oscillator con- 
figuration 19] is based on the two-resonance picture for 
weakly non-integrable systems The strongest 

resonant term in the non-integrable perturbation of an 
integrable system (the so-called guiding resonance) pro- 
duces a separatrix in the system phase space. This sepa- 
ratrix is destroyed by another term in the perturbation, 
which creates a thin stochastic layer in the surrounding 
part of the phase space. In a disordered system, the 
main contribution to the chaotic phase space comes from 
those configurations of disorder and from those regions of 
the phase space where both perturbation terms are res- 
onant, so one cannot really distinguish between the one 
responsible for the appearance of the separatrix and the 
one responsible for its destruction. It is crucial that two 
resonance conditions should be met simultaneously. 

Resonances involving many oscillators are expected to 
give a subleading contribution to the chaotic phase vol- 
ume, as the corresponding perturbation terms can be 
generated in high orders of perturbation theory, thus re- 
sulting in high powers of e. The minimal number of os- 
cillators needed to generate chaos in the model of Eq. ([1]) 
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is two [24j |. However, even if u]% ~ ll>2, the frequency 
of the separatrix-destroying perturbation, wi + u> 2 , can- 
not be small for the chosen disorder distribution, 1/2 < 
ijj 2 n < 3/2, which leads to an exponential suppression of 
the chaotic phase volume (l#jlo [. In fact, for the par- 
ticular case of Eq. ([1]), the situation is even worse: no 
separatrix exists in the phase space of the slow motion 
(see Appendix I A II) . Thus, the dominant contribution 
to the chaotic phase volume should come from triples of 
oscillators with all three frequencies close to each other. 
The oscillators should also be neighboring each other in 
space, since coupling distant oscillators requires high or- 
ders of perturbation theory and results in a high power 
of e. Thus, w(e) is essentially determined by the prob- 
ability (per unit length along the chain) to find three 
neighboring oscillators whose frequencies differ by ~ e. 
This fixes the power w(e) oc e 2 . 

To put this argument on a quantitative basis, we as- 
sume u)\ s» u>2 « 0J3, and average Eq. (JTJ) over fast os- 
cillations (see Appendix IA 21 for details). This gives the 
effective Hamiltonian of the resonant triple: 

H tr = r!|*!| 2 - (n + n')\v 2 \ 2 + o'|* 3 | 2 + 

+ \ (l*il 4 + 1* 3 | 4 + l*i - * 2 | 4 + 1*2 - * 3 | 4 ) , 

(7) 

written in terms of complex canonical variables i^* , 
related to p n , q n as 



(8) 

where u) = (uj\ +u} 2 + W3) /3. The rescaling of \l/'s by y/1 
restricts them to the unit sphere, |*i| 2 +|*2| 2 +|*3| 2 = 1, 
which is invariant under the dynamics generated by the 
Hamiltonian in Eq. (JTJ)- The rescaled frequency mis- 
matches 0, ft' are defined as 



^ _ 2uji -u) 2 -u} 3 ^, 2u 3 —W2-U1 



(9/4)(I/a>») 



(9/4)(//u, 2 ) ' 



(9) 



For the fixed-e ensemble, Eq. (j4)), we have / = 
3e/cl> + 0(e 2 ), and the initial condition should be cho- 
sen in the form *i = #2 = e^+^'/v^ 
\&3 = e~ lv /\/3, where < tp, ip' < 2n are random phases 
(the global phase drops out of the result). Then, defining 
6(<p, f'; n, ft') to be 1 if the corresponding trajectory is 
chaotic for given values of ip, ip' ,fl,fi' and otherwise, 
we obtain from Eq. 

OO 27T 

w(e) = 8u 2 { dfidn' [ ^e(yy;fi,no. (10) 

J J 2tt 2n 



Since neither the Hamiltonian [Eq. ([7])] nor the initial 
condition contain any small parameters, the integral 
over ft, ft' in Eq. flTDJl is dominated by |fi'| ~ 1, 
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FIG. 1. A grayscale plot of the integral 

J^ w {dif/2n)(d(p' /2tt) Q{ip, ip'; 0, Q!) as a function of Q, ft'. 



while for larger mismatches the chaotic regions quickly 
shrink. Thus, the limits of the ft, ^'-integration (which 
are ~ 1/e) have been extented to infinity. The factor e 2 in 
front of the integral appears because fi, Of oc 1/1 oc 1/e, 
see Eq. ©. The integral is evaluated numerically to be 
16.9±0.2, leading to w(e -> 0) = Ae 2 with A = 1.37- 10 3 , 
as stated in Eq. ([6]). 

The key element of the numerical procedure is the cri- 
terion which enables one to distinguish between regular 
and chaotic motion. Formally, one studies the largest 
Lyapunov exponent er, characterizing the mean exponen- 
tial rate of divergence of two initially close trajectories, 



lim 



A(t) 



t 



A(t) = In 



d(oy 



(11) 



where d(t) is the distance between points belonging to 
the two trajectories [l6j]. If a = 0, the motion is regular, 
if a > 0, it is chaotic. In practice, one may integrate the 
equations of motion for a sufficiently long time T, and 
consider a = if A(T)/T < (a few times)(l/T) dflU^. 
Here we use a slightly different criterion, based on the 
behavior of A(f) in the whole integration interval < t < 
T. Namely, we check how well A(t) can be approximated 
by a logarithmic function, as described in Appendix [B] 

In Fig. [TJ the <p, (//-integral in Eq. (fTU|) is plotted as a 
function of ft, ft' . The fact that the chaotic region is con- 
fined in all directions, represents a numerical proof of the 
statement that chaos arises mostly in the regions where 
two resonant conditions are satisfied simultaneously. 

Using the same algorithm, w(e) was also calculated 
by direct numerical sampling of disorder realizations on 
chains of lengths L = 8,16,32,64,128. For e > 1CT 4 , 
the integration time T = 10 6 and averaging over 10 4 
realizations was performed for each L; for smaller e < 
10~ 4 , when the probability Pz,(e) is also very small, more 
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FIG. 2. (color online) (1/L)|ln(l - Ph{e))\ versus e for 
the fixed-density excitation and different values of L = 
8,16,32,64,128 (symbols). The error bars for L = 128 cor- 
respond to the relative error l/y/N, where N is the absolute 
number of detected events. The dashed line represents the 
dependence w(e) = (1.37- 10 3 ) e 2 , obtained from the resonant- 
triple calculation. 
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FIG. 3. (color online) The probability of chaos versus e for a 
single-site excitation and two values of L — 11,21 (symbols). 
The dashed line represents the dependence 235 e 2 , obtained 
from the resonant-triple calculation of Appendix [C] 



IV. RELATION TO WAVE PACKET 
SPREADING 



than 4 • 10 4 realizations were required to accumulate at 
least N — 20 events for each data point. Then 1/yN was 
assumed to give the relative uncertainty of the obtained 
value w(e). Also, for e < 10~ 4 , the integration time 
had to be increased to T = 3 ■ 10 6 in order to reliably 
distinguish between regular and chaotic dynamics (see 
also Appendix |B|) . 



The numerical results, as well as the dependence 
w(e) = (1.37 • 10 3 ) e 2 , obtained from Eq. (HO), are shown 
in Fig. [2] by the symbols and the dashed line, respec- 
tively. Fig. [2] represents the main result of the present 
work. Starting from L = 16, the collapse of the numeri- 
cal data is very good, which numerically proves the local 
origin of chaos. Ate<3-10~ 4 the numerical data fall 
on the dashed line. At larger e the data collapses on 
some intermediate asymptotics for w(e), which can be 
approximated by a power law w(e) — Be 13 with the ex- 
ponent P = 2.85 ± 0.1 and a surprisingly large prefactor 
B = (1 - 5) • 10 6 . The precise mechanism responsible 
for this intermediate asymptotics is not clear at the mo- 
ment. This mechanism is likely to involve a larger num- 
ber of oscillators (more than three), since (i) it is in this 
region of e that the symbols for short chains (L = 8) ex- 
hibit a systematic deviation from those for long chains in 
Fig. [2] in the direction of suppression of the intermediate 
asymptotics; (ii) for a single-site initial excitation (see 
Fig. |3] below) the intermediate asymptotics is completely 
absent, indicating that it requires a certain initial spread. 



In this context, one should mention the result reported 
in Ref. !22] , where the probability of subdiffusive spread- 
ing of an initially localized wave packet was studied for 
exactly the same model, Eq. ((T). Namely, it was argued 
that an initially localized wave packet either stays local- 
ized forever or spreads indefinitely, and the probability 
of spreading was estimated to be proportional to e at 
small e. In the present work, the probability of chaos, 
was shown to scale as e 2 . Since the latter is smaller 
than the former, the combination of the two results would 
suggest that the wave packet spreading does not require 
chaos, which is quite counter-intuitive. However, in the 
author's opinion, it is more plausible that the probability 
of spreading was strongly overestimated in Ref. [2^, as 
discussed below. 

First, let us exclude the trivial possibility that the dif- 
ference between the two results is simply due to the fact 
that initial conditions studied in Ref. [22} correspond to 
excitation localized on a single site, while a finite density 
was assumed in the present work. Indeed, for a finite 
density the trajectory is chaotic when a resonant triple 
occurs anywhere on the chain, while for a single-site ex- 
citation the condition is simply that the resonant triple 
includes the excited site. This affects the scaling with L, 
but not with e: the probability of chaos for a single- 
site excitation scales as Pi(e) = A\t 2 and is independent 
of L. The resonant-triple calculation gives A\ w 235 
(see Appendix [C), which is also confirmed by the direct 
numerical sampling, as shown in Fig. [3] (Note also the 
direct crossover from 1 to A\e 2 without any intermediate 
asymptotics). 

The analytical arguments of Ref. [22[ are based on 
the assumption that a single resonance is sufficient for 
spreading. This immediately gives the probability seal- 



5 



0.2 - 



0.15 



0.1 - 



0.05 




FIG. 4. (color online) Partial contributions to the probability 
distribution of Q = 1/(V — 1) [see Eq. (|12p ] from regular and 
chaotic trajectories (upper and lower curves, respectively) at 
times t = 10 5 , 10 6 [thick red (gray) and thin black curves, 
respectively] for two chain lengths L — 11,21 (dotted and 
solid curves, respectively) for e = 0.01. 



ing oc e. However, it is not clear why a single resonance 
should lead to unlimited spreading; indeed, an isolated 
nonlinear resonance is known to produce just periodic 
oscillations 0, EH- 

The numerical procedure of Ref. [22[ was based on the 
study of the participation number, 



(12) 



where the on-site energy corresponding to the model of 
Eq. fl} can be defined as 



p _ pI 

&n - TT- 



+ o [0?n - <7n-i) 4 + (q n - q n +i) 4 



for all n except n = 1 and n = L, where only one nonlin- 
ear term corresponding to the unique neighbor should be 
taken. For a perfectly thermalized chain V = L/2, and 
for a state perfectly localized on a single site V = 1. In 
Ref. [13 , a trajectory was counted as spreading if V ex- 
ceeded an arbitrarily chosen threshold of 1.2 by the time 
t = 10 9 . Thus, the main assumption behind the numerics 
was that if a trajectory has overcome the limit V = 1.2 
at t = 10 9 , it will spread forever. In the following, it is 
argued that this is unlikely to be true for the majority of 
the trajectories. 

Let us analyze the probability distribution of Q = 
l/(V — 1), which is a more convenient quantity to an- 
alyze at small e, when the spreading trajectories corre- 
spond to Q ~ 1, while the overwhelming majority of 
strongly localized solutions form the broad large- Q tail. 
In Fig. @] we separately plot the two contributions to the 
distribution of Q from regular and chaotic trajectories 
at t = 10 5 , 10 6 for L = 11.21. The value of e = 0.01 is 



sufficiently small to be in the asymptotic regime, as seen 
from Fig. [3J and more than 50000 disorder realizations 
were used for each of the four pairs of curves. The differ- 
ence between L = 11 and L = 21 curves is not detectable 
within the numerical resolution, while a slight difference 
between t = 10 5 and t = 10 6 curves can be seen on the 
low-Q side of the chaotic curves, which indicates some 
spreading of the wave packet. The fact that a large part 
of the probability distribution at Q < 5 (corresponding 
to V > 1.2) belongs to the large-Q tail due to localized 
trajectories, suggests that many trajectories counted as 
spreading in Ref. (22[, in fact were not. 

Another argument, not relying on numerical integra- 
tion over long times, can be given by considering a regular 
solution, predominantly localized on site n. The ampli- 
tude on the neighboring sites n — 1, n + 1 can be found 
by perturbation theory, 



q n = cosw„t + 0(e 3/2 ), 



e 3 / 2 / 3cosu; n £ cos3u;„£ 



9n±l 



\ w n±l - ^1 W 2 ±1 - Qui 



(14a) 
(14b) 



which gives the participation number 



V w 1 



9e 2 \^ 



w n±l 



J n± 



i/9 



K ±1 - ul) 2 K 



±i 



9^ 2 



(15) 



Here we omitted oscillating terms, which at large t 1 
quickly vanish upon disorder averaging), as well as terms 
which do not contain potentially small denominators and 
thus are bounded by O(e) [such as 3e/(4o;^), originating 
from the (q n ±i — qn) 4 contribution to E n ±i\. 

From Eq. dl~5l) . one can find the probability of V > 
1.2 at small e by noting that V > 1.2 may occur when 
either w„+i or uj n -\ is close to u) n . The probability of 
both occurring simultaneously can be neglected, as well 
as the probability of uj n ±i being close to 3w ra , as these 
probabilities are oc e 2 . Then, the probability of V > 1.2 
is given by 



(w 2 ±1 - c4)2 



2Pr<^ \ojI 



J n+1 



> 0.2 



3e 



/02 s/2ul 



3 /2 

3\/8e f duil 




20.8 e. 



(16) 

The coefficient in front of e is very close to that obtained 
numerically for the probability of spreading (see inset of 
Fig. 3 of Ref. [13). 
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The fact that Eq. flTSl) . derived under the assumption 
of validity of perturbation theory at arbitrarily long times 
(i. e., that the trajectory is localized), gives practically 
the same result as the numerics of Ref . [22| , means that 
almost all trajectories counted as spreading in Ref. [22j . 
were in fact localized. Thus, the question of whether the 
probability of unlimited wave packet spreading scales as 
oc e, or oc e 2 , or is identically zero, remains open. 



V. CONCLUSIONS 

In this work, the probability for a long chain of os- 
cillators, defined by its Hamiltonian, Eq. (fTJ), to be on 
a chaotic trajectory, was calculated by analyzing the 
chaotic phase space of rare resonant configurations of 
three oscillators. The result, 1 — e~ At L , agrees with the 
direct numerical evaluation of Lyapunov exponents for 
many disorder realizations and initial conditions at suffi- 
ciently small values of e, the energy per oscillator. While 
the coefficient A depends on the specific ensemble used 
for the initial conditions, the power law e 2 is determined 
by the mere fact that chaotic phase volume is dominated 
by the regions of the phase space where two resonant 
conditions for three neighboring oscillators are satisfied 
simultaneously. This fact was numerically verified for the 
three-oscillator configuration, see Fig. [TJ 

The following main conclusions can be drawn from the 
results of present work, (i) The local nature of chaos, 
established earlier for the discrete nonlinear Schrodinger 
equation with disorder (DNLS) in Refs. pjl [2(| for weak 
nonlinearity, is further confirmed here for a different 
model. Namely, the probability of regular dynamics de- 



cays with length L as e~ wL , so w can be called the prob- 
ability per unit length for the chaos to occur, (ii) Two 
different ways to calculate this probability, namely, from 
the analysis of the phase space of a resonant triple, and 
by direct numerical sampling of disorder realizations, give 
the same result, but this agreement sets in at unexpect- 
edly small values of e ~ 3 x ICF 4 . At larger e, an interme- 
diate asymptotics is seen in the numerical results. This 
fact may be relevant for understanding the disagreement 
between the results of Refs. [l9| and [2(| for the chaotic 
probability in the DNLS; indeed, one cannot exclude exis- 
tence of similar intermediate regimes in DNLS. However, 
the results of the present work were obtained for a quite 
different model with a different number of independent 
parameters, so no quantitative comparison can be made 
with the DNLS, and the latter should be studied further. 
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Appendix A: Few-oscillator configurations 

It is convenient to pass to the canonical action-angle 
variables, (p n ,q n ) -> {I n ,4>n)'- 



p n = y/2uj n I n cos0„, q n = \ — sin^„. (Al) 

V w n 

Hamiltonian (TT]), averaged over the fast oscillations, (that 
is, with terms depending on (j) n + (j) n+ i omitted) takes the 
following form: 



L-l 



H 



In , In+1 
W n <^n+l 



/ I n I n +l , 

2t cos( 

V w n uj n +i 



n +l) 



(A2) 



It should be noted that the part of the Hamiltonian cor- 
responding to two or three oscillators (say, n, n+1, n+2) 
far from the ends of the chain is different from the Hamil- 
tonian given by Eq. (|A2|) for L = 2 or L = 3. Indeed, 
in the former case the Hamiltonian of the nth oscillator 
contains the term (3/4)/ 2 /w 2 , coming from coupling to 
the two neighbors n—1 and n+1, while in the latter 
case the n = 1 oscillator has only one neighbor, so the 
nonlinear term is twice smaller, (3/8)/ 2 /u; 2 . 



Two oscillators 



Let uj n m uJn+i be strongly different from u) n -i, to n +2- 
We perform a canonical change of variables of the pair 



I — In + In+1, 
j In — In+1 



L>n+1 



(A3) 



J n+1, 



7 



and further denote J = (7/2) cos??. Then the Hamilto- 
nian of the pair is given by 



whose conjugate momenta are itp*, iip2, *^3i respectively, 
we can write the Hamiltonian of the triple as 



H = "1+^2 / + ( w i - cj 2 )(//2) cost) + 



3 I ( UJ\ + U>2 U>1 — UJ2 



4 4 \ U!lU>2 

3 (7 2 /4)sintf 

2 ( W iw 2 )3/ 2 

3 I 2 / 4 



COS I? ) + 



U>lOJ 2 

[(jJl + 0J2 — (0Jl — W 2 ) COS 1?] COS y> + 



4 W1W2 



sin 1? cos 2ip, 



(A4) 



where = w„ + 4J„_i/(w„a; n _i), w 2 = + 
4/ n+2 /(u; n+1 u; n+2 ), We can assume i n _i, / n + 2 to be con- 
stant, as their changes due to weak non- resonant pertur- 
bation are small. This Hamiltonian can be approximately 
rewritten as 

3 I 2 

H fa — — (f2 cos d + 4 sin •& cos <y9 + sin 2 $ cos 2w) + 
16 uj l 



const, 
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■ w 2 



(A5a) 
(A5b) 



This Hamiltonian has only two elliptic stationary points 
in the phase space (1?, <p) at any value of the rescaled 
mismatch fi. Hence, its phase space does not contain a 
separatrix at all, so we do not expect any chaos at small e. 



2. Three oscillators 

As we are interested in the properties of a resonant 
triple n, n + 1 , n + 2 far from the ends of the chain, let us 
denote the three frequencies by 

4/n-l 

w x = uj n H , 



U>2 — Wn+1) 
/ 

W 3 = OJn+2 + 



(A6) 



47, 



n+3 



and assume I n -i and in +3 to be constant, since their 
variation is small in the parameter e -C 1. Moreover, 
we can also neglect the modification in the probability 
distributions of uj[,uj 3 with respect to that of u) n , as the 
difference is again ~ e. The most important contribution 
to the chaotic fraction will come from the region where 
the differences uj[ — uj' 2 , w 2 — uj 3 , are small compared to 
the frequencies themselves. Having this in mind, let us 
further denote 

u = -1 i K (A7) 



Introducing the complex canonical coordinates ip\ — 

~i<t>n ^ 2 = ./TTT-, p-i<t>n+l ,/,, = ./L,, p-i<t>n + 2 



3 i>*i>*ipiipi 3 ^* 3 ipiii 3 ^ 3 



K) 2 8 K) 2 



8 2^ 



1 vX^. 



1 



(A8) 



In the nonlinear terms (those of the fourth order in 
tpi, ip*) it is sufficient to replace 2 3 —> w, to the leading 
order in e. The total norm, 



i = ^*V>1 + ^2^2 + ^3 ^3, 



(A9) 



is conserved for the Hamiltonian (|A8I) . and the corre- 
sponding conjugate phase, 



(A10) 



even though depends on time in some complicated way, 
does not contribute to anything. In the fixed- e ensemble, 
we simply have I = 3e/ui + 0(e 2 ). It is convenient to pass 
to rescaled variables 



Vi = -jjil>i, i = 1,2,3, 

4 CO 2 2i0\ — L0' — LOn 

n = — - — -, 

3 7 3 

, 4 Co 2 2a4 — cjo — u\ 

Q = 2 2 i-. 

3 I 3 



(All) 



Then |*i| 2 + |* 2 | 2 + |* 3 | 2 = 1. The probability density 
for frequencies is 



3/2 

v(ui) = lim / dxdy dz 5 \ - 
x 1 n,«'->o J \ 



1/2 



(9/4)(i/^ 2 

' (9/4)(//^ 2 ) 
+ 0(e 2 ). 



27 P 

y 7d 



(A12) 



The rescaled Hamiltonian is given by 



(3/4)i 2 c 



n^xl 2 -(n + Q,')\v 2 \ 2 + n'\y 3 \ 2 



1 



1 



2 1*11 +^1*31 + 

i|*l-* 2 | 4 + i|* 2 -4'3| 4 . 



(A13) 
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For the thermal distribution the initial condition should 
be taken in the form 



*2 = sjx-ye 



^3 = VV 



(A14) 



Then the chaotic fraction is given by 



w(e) = I i 
o 



-ul/e 



Pdl 



27 I 2 

did I dfl dO' x 

2 Co 



'1/2 



X / ~2n~2^ J dX J dy0 ( a; ' y;i/3l '^ 2; ^'^')- 



(A15) 

The first two integrals amount to 216 e 2 . The fixed-e 
ensemble corresponds to fixed x = 2/3, y — 1/3, so the 
chaotic fraction is given by 




J 2 



J, 




FIG. 5. The allowed region of Ji, J2, shown by the shaded 
area. 



'3/2 



w(e) 



27 (3e) 



2 Co 3 



duo I dSldfl' x 



1/2 

2tt 



(A16) 



/ ©(VS. 1/3; 



The first integral is equal to 81 e 2 , and the integral over 
fl,Q' to 16.9 ±0.2. 

One can separate the conserved total norm and write 
the explicit Hamiltonian for the two remaining degrees of 
freedom. Let us label the actions and phases of the oscil- 
lators n, n + 1, n + 2 by i = 1,2,3, respectively. Then the 
corresponding canonical transformation is determined by 



+ 92 + 93 



I T 2I x -I 2 -h 

2^1 = ^ > V\=9\ 

I h + h~2h 

gJ2 = g , f2 = (P2 



(A17) 



To ensure I\, l2,h > 0, the variables J\, J2 should lie 
inside the triangle, shown in Fig. [5l The inverse trans- 
formation is 



1p! = y/(l + Ji)//3 e < (- 2 *'i+V2)/3-i* ) 

^2 = y/(l -Ji- J2)//3e 4(vi+¥,2)/3 -^, (A18) 

1P3 = + J2)I/3e l{lpl ' 2<P2)/3 ' 1 ' 1 '. 

The rescaled Hamiltonian for the degrees of freedom 



(J 1 ,ip 1 ,J 2 ,^2) is given by 

— — - = (20 + fi' - 2/3) Ji + (0 + 20' - 2/3) J 2 
(//2cj) 2 



v/(l + Ji)(l- Ji- J 2 )(2- J 2 )cos 



+ q i 1 + J i)( 1 ~ J i ~ J2)cos2v?i 



J 2 )(l- Ji+ J 2 )(2- Ji)cos 



¥>2 



+ o (1 - ^2)(1 - Jl + J2) COS2(^ 2 

O 

7 2 

+ 3-3 JlJ - 



(A19) 



Still, due to the presence of trigonometric functions, the 
numerical integration of the equations of motion for this 
Hamiltonian would be less efficient than for the polyno- 
mial Hamiltonian ([7J with three degrees of freedom. 



Appendix B: Numerical criterion for chaos 



The numerical criterion for the chaotic motion used in 
the present work is based on the analysis of A(t) = In , 

where d(t) is the distance between points belonging to 
two initially close trajectories. To illustrate the idea, 
we plot several traces A(t) in Fig. [5] Looking at them, 
it is quite easy to guess which curve corresponds to a 
regular motion, and which one to chaotic. The only ex- 
ception is the lowest curve in Fig. HJd). We know it 
should be chaotic, as the corresponding realization of 
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D[A(t)] 
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time (10 6 ) 
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0.2 0.4 0.6 0.8 1 

time (10 6 ) 




FIG. 6. (color online) A(t) for < t < 10 6 . Panels (a), (b), (c) correspond to three chain lengths: (a) L = 3, (b) L — 10, 
(c) L — 32. All curves on a given panel correspond to the same realization of the disorder {w^}, phases {4>n}, and differ only 
by the value of e n — e. Higher curves correspond to larger values of e. Panel (d) corresponds to L = 32, but the disorder 
realization was artificially modified by setting UJ12 = W13 = W14, thereby creating two neighboring resonances. 



{w^} contains two neighboring resonances, introduced 
"by hands". However, at e = 10 ~ 5 energy exchange 
between the oscillators is so slow, that following their 
dynamics up to T — 10 6 is insufficient to distinguish 
between regular and chaotic motion. This gives us the 
lower boundary on e. Compare also the lowest curve in 
Fig. [He) and the middle curve in Fig. [HJd). The former 
one is regular and the latter is chaotic, even though they 
have close values of A(t = T). Thus, the fixed cutoff for 
the Lyapunov exponent (independent of L and e) , which 
was used in Ref. [2(|, as the criterion for chaos, may in- 
troduce systematic errors and affect the scaling. 

The main observation from Fig. [6] is that for chaotic 
trajectories A{t) grows linearly on the average, corre- 
sponding to a finite Lyapunov exponent, while for regular 
trajectories 

A(t) = \nt + q , (Bl) 

up to some weak noise, qo being some realization- 
dependent constant. This can be understood very sim- 



ply: if two trajectories lie on neighboring invariant tori, 
their frequencies are only slightly different, so the phase 
mismatch between them is accumulated slowly and lin- 
early in time. The distance between the trajectories is 
proportional to this phase mismatch, as long as the latter 
is small compared to unity. 

To determine how well a given function A(t) is ap- 
proximated by the form (jBll) on an interval < t < 
T, consider the time average (denoted by overline): 
[A(t) — hit — q] 2 . It is a quadratic function of the pa- 
rameter q, which has a minimum at some q, and its value 
at the minimum is D > 0. If A(i) is given exactly by 
Eq. (|Bll . then simply [A(t) - hit - q} 2 = (q - q ) 2 , so 
D = 0. Thus, the value at the minimum, given by 

D[A(t)\ = \Mj~) - Int} 2 - (JA(<) - ln£]) 2 , (B2) 

measures the quality of the fit of A(i) by the expres- 
sion (jBll) . D[A(t)] for all curves in Fig. [5] is shown on 



the corresponding panel. Thus, we consider the motion 
as regular if D[A(t)} < 1 and chaotic otherwise. 

When the Lyapunov exponent a > 0, it sets the natu- 
ral upper time limit for the numerical integration of the 
differential equations. Indeed, when the double-precision 
machine zero, 2~ 52 , multiplied by e CTt becomes of the or- 
der of unity, the integration inevitably deviates from the 
original trajectory, no matter how good the integration 
scheme is. This corresponds to A = ln2 D2 ss 36. To 
check this, we have performed the time-reversal test: at 
time t = t*, we invert all the momenta, and integrate 
up to t = 2t*. Ideally, the system should return to 
its initial condition. This was found to be the case if 
A(i*) < 25 — 30 (for many trajectories whose Lyapunov 
exponents differed by one or two orders of magnitude). 
Thus, considering larger A would produce the Lyapunov 
exponent not for a given trajectory, but its certain av- 
erage over the phase space. As we are interested not in 
the trajectory itself, but only in whether it is chaotic or 
regular, we consider the very fact that A(t) has reached 
25 a sufficient evidence for chaos. If this happens within 
the integration time, t < T, we stop the integration, and 
simply set D[A] = 30. 

One could think that for every realization of the disor- 
der {w 2 } and of the oscillator phases {<fi n } exists a thresh- 
old value e c , such that for e < e c the motion is quasiperi- 
odic, and for e > e c it is chaotic. Indeed, a natural guess 
is that upon increase of the integrability-breaking param- 
eter e the chaotic region of the phase space should grow. 
This would make the calculation for the fixed-e ensemble 
more efficient, as one would not need to probe too small 
or too large values of e. However, this guess turns out to 
be wrong, as seen from Fig. [7J where we plot lnD[A] as 
a function of e in the vicinity of e = 10~ 3 for the same 
realization of disorder and phases as in Fig. life). The 
flat regions on the level of 3.4 ... = In 30 correspond to 
A(t) reaching 25 for t < T, as discussed in the previous 
paragraph. One can see that the conclusion about reg- 
ular/chaotic character of the dynamics is little sensitive 
to the chosen numerical border D = 1: had we chosen 
D = 5 or D = 1/5, the result would not change. Wc 
also note that for the particular realization, correspond- 
ing to Fig. [JJ in all intervals of e where the dynamics 
is chaotic, the Lyapunov eigenvector is confined to sites 
from n = 21 to n = 26. This is clearly related to the fact 
that in this realization \uj22 — I w 0.9 • 10 -3 , so the 
observed reentrant behavior of chaos seems to occur for 
the same guiding resonance. 

To check how the results are affected by the choice of 
the integration time T, we take the same realization of 
disorder as for the trajectories in Fig. [5] (the chain length 
L = 32), and set "by hands" three nearby frequencies to 
be 

9 e 9 e 

U3i2 = uj+- — fi, u 13 = u-- — Q, u)h = oj, (B3) 
4 w J 4 w J 

with a fixed <D = v(L7, two values of e = 10~ 4 , 0.3 ■ 10~ 4 , 
and O varying in the interval — 2 < < 6. This chain 



< 
c 




6(10 ) 

FIG. 7. In 73 [A] [with D defined in Eq. <[E2)l] as a function of e 
for the same realization of disorder and phases as in Fig. Oc) . 
The upper cutoff corresponds to In 30, as discussed in the text. 



is supposed to be well described by the effective Hamil- 
tonian in Eq. ([JJ with Q! = 0. On the two panels of 
Fig. [5] we plot the phase-averaged chaotic fraction as a 
function of fi for the two above-mentioned values of e, 
and for T = 10 6 , 3 • 10 6 , 5 • 10 6 , together with the re- 
sult of the calculation using Eq. ([JJ (the section of Fig. Q] 
along the line fi' = 0). For both values of e, as T is in- 
creased, the curves approach the result of Eq. ([JJ, up to 
an overall horizontal shift. This horizontal shift appears 
because Eq. (|B3j) neglects the nonlinear frequency shift 
for the oscillators with n = 12 and 14 due to their cou- 
pling to the oscillators with n — 11 and 15, respectively 
[the difference between primed and unprimed frequen- 
cies in Eq. (|A6|) ]. This is also the reason why the dip 
at = for the result of Eq. ([JJ is not resolved on the 
curves for the long chain. What is important, is that 
while at e = 10~ 4 the integration time T = 10 6 is suffi- 
cient, for e = 0.3 ■ 10~ 4 integration up to T = 10 6 clearly 
underestimates the chaotic fraction, so that T = 3 • 10 6 
is required. 

Finally, we mention that the numerical integration of 
the differential equations was performed using (i) the 
fourth-order Runge-Kutta method for the calculations 
reported in this section, and (ii) Bulirsch-Stoer method 
with polynomial extrapolation for the accumulation of 
statistics. The latter method turns out to require about 
ten times less CPU time that the former to reach the 
same level of accuracy. The energy conservation was sat- 
isfied extremely well in all calculated trajectories, includ- 
ing those with largest Lyapunov exponents, even at times 
when the time-reversal test had failed. Thus, it is unlikely 
that use of an integration scheme which automatically re- 
spects some conservation laws of the original equations 
(e. g., a symplectic integrator which preserves the phase 
space volume) would lead to more accurate results for 
a given trajectory: the accuracy is lost primarily in the 
directions, orthogonal to conservation laws. 
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Appendix C: Single-site excitation 

Let us determine the probability of chaos for an initial 
condition, localized on a single site, 

p„(0) = 0, q n (0) = 5 n} ( L+ iy 2 V2e/u n , (CI) 

for L odd. Chaos should occur if the resonant triple 
contains the site n = (L + l)/2. Thus, the probability 
of chaos Pi(e,L) should not depend on L and scale as 
Pi(e) = Aie 2 , at large L and small e. 

To determine the coefficient A\, we again use the ef- 
fective Hamiltonian of the triple, Eq. (JT]), and consider 
two initial conditions: 



*i(0) = -^=, * 2 (0) - * 3 (0) = o 



*i(0)-0, * 2 (0) = 



V3' 



(C2a) 
* 3 (0) = 0. (C2b) 



i. e., the excitation to be initially localized on one of 
the lateral sites of the triple, or on the central site. Let 
9i i2 (^,ri') equal one if the corresponding trajectory is 
chaotic and zero if it is regular, for each of the two initial 
conditions, respectively. The functions 0i,2(^,f^') are 
shown in Fig. [3] The chaotic probability is then given by 



FIG. 8. (color online) The phase-averaged chaotic fraction 
P(fi) for a given realization of disorder on a chain with L = 32 
and frequencies U12, ujis, CJ14 chosen according to Eq. (|B3I) . 
for (a) e = 1(T 4 , (b) e = 0.3 ■ 10~ 4 . The thick solid line on 
both panels corresponds to the calculation using the effective 
Hamiltonian from Eq. J7J and represents the section of Fig. Q] 
along the line Q' = 0. The thin solid, dashed and dotted lines 
correspond to T = 10 6 , 3 • 10 6 , 5 ■ 10 6 , respectively [the last 
curve is present only on panel (b)]. 



oo 

P 1 (e) = 81e 2 y dOdn' [2 01(0,0') + 9a(n,n')]- ( C3 ) 



For the initial condition (|C2al) . the ft, O-integral is equal 
to 0.75±0.02, while for Eq. (jC2b|) it is equal to 1.39±0.02, 
which gives A\ sa 235. 
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